home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / min_curve_surf.pro < prev    next >
Text File  |  1997-07-08  |  10KB  |  277 lines
  1. ; $Id: min_curve_surf.pro,v 1.8 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1993-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5.  
  6. FUNCTION min_curve_surf, z, x, y, REGULAR = regular, XGRID=xgrid, $
  7.     XVALUES = xvalues, YGRID = ygrid, YVALUES = yvalues, $
  8.     GS = gs, BOUNDS = bounds, NX = nx0, NY = ny0, XOUT=xout, $
  9.     YOUT=yout, XPOUT=xpout, YPOUT=ypout, TPS=tps
  10. ;
  11. ; Copyright (c) 1993, Research Systems, Inc.  All rights reserved.
  12. ;    Unauthorized reproduction prohibited.
  13. ;+
  14. ; NAME:
  15. ;    MIN_CURVE_SURF
  16. ;
  17. ; PURPOSE:
  18. ;    This function Interpolates a regularly- or irregularly-gridded
  19. ;    set of points with either a minimum curvature surface or
  20. ;    a thin-plate-spline surface.
  21. ;
  22. ; CATEGORY:
  23. ;    Mathematical functions, Interpolation, Surface Fitting
  24. ;
  25. ; CALLING SEQUENCE:
  26. ;    Result = MIN_CURVE_SURF(Z [, X, Y])
  27. ;
  28. ; INPUTS: 
  29. ;    X, Y, Z:  arrays containing the X, Y, and Z coordinates of the 
  30. ;          data points on the surface. Points need not be 
  31. ;          regularly gridded. For regularly gridded input data, 
  32. ;          X and Y are not used: the grid spacing is specified 
  33. ;          via the XGRID and YGRID (or XVALUES and YVALUES) 
  34. ;          keywords, and Z must be a two dimensional array. 
  35. ;          For irregular grids, all three parameters must be
  36. ;          present and have the same number of elements. 
  37. ;
  38. ; KEYWORD PARAMETERS:
  39. ;    TPS:    If set, use the thin-plate-spline method, otherwise
  40. ;        use the minimum curvature surface method.
  41. ;  Input grid description:
  42. ;    REGULAR:  if set, the Z parameter is a two dimensional array 
  43. ;          of dimensions (N,M), containing measurements over a 
  44. ;          regular grid. If any of XGRID, YGRID, XVALUES, YVALUES 
  45. ;          are specified, REGULAR is implied. REGULAR is also 
  46. ;          implied if there is only one parameter, Z. If REGULAR is 
  47. ;          set, and no grid (_VALUE or _GRID) specifications are 
  48. ;          present, the respective grid is set to (0, 1, 2, ...). 
  49. ;    XGRID:    contains a two element array, [xstart, xspacing], 
  50. ;          defining the input grid in the X direction. Do not
  51. ;          specify both XGRID and XVALUES. 
  52. ;    XVALUES:  if present, XVALUES(i) contains the X location 
  53. ;          of Z(i,j). XVALUES must be dimensioned with N elements. 
  54. ;    YGRID:    contains a two element array, [ystart, yspacing], 
  55. ;          defining the input grid in the Y direction. Do not
  56. ;          specify both YGRID and YVALUES. 
  57. ;    YVALUES:  if present, YVALUES(i) contains the Y location 
  58. ;          of Z(i,j). YVALUES must be dimensioned with N elements. 
  59. ;
  60. ;  Output grid description:
  61. ;    GS:      If present, GS must be a two-element vector [XS, YS],
  62. ;          where XS is the horizontal spacing between grid points
  63. ;          and YS is the vertical spacing. The default is based on
  64. ;          the extents of X and Y. If the grid starts at X value
  65. ;          Xmin and ends at Xmax, then the default horizontal
  66. ;          spacing is (Xmax - Xmin)/(NX-1). YS is computed in the
  67. ;          same way. The default grid size, if neither NX or NY
  68. ;          are specified, is 26 by 26. 
  69. ;    BOUNDS:   If present, BOUNDS must be a four element array containing
  70. ;          the grid limits in X and Y of the output grid:
  71. ;          [Xmin, Ymin, Xmax, Ymax]. If not specified, the grid
  72. ;          limits are set to the extent of X and Y. 
  73. ;    NX:       The output grid size in the X direction. NX need not
  74. ;            be specified if the size can be inferred from GS and
  75. ;          BOUNDS. The default value is 26.
  76. ;    NY:       The output grid size in the Y direction. See NX. 
  77. ;    XOUT:      If present, XOUT must be a vector containing the output
  78. ;          grid X values. If this parameter is supplied, GS, BOUNDS,
  79. ;          and NX are ignored for the X output grid. Use this
  80. ;          parameter to specify irregular spaced output grids.
  81. ;    YOUT:      If present, YOUT must be a vector containing the output
  82. ;          grid Y values. If this parameter is supplied, GS, BOUNDS,
  83. ;          and NY are ignored for the Y output grid. Use this
  84. ;          parameter to specify irregular spaced output grids.
  85. ;    XPOUT:      If present, the arrays XPOUT and YPOUT contain the X and Y
  86. ;    YPOUT:      values for the output points. With these keywords, the
  87. ;          output grid need not be regular, and all other output
  88. ;          grid parameters are ignored. XPOUT and YPOUT must have
  89. ;          the same number of points, which is also the number of
  90. ;          points returned in the result.
  91. ;
  92. ; OUTPUTS:
  93. ;    This function returns a two-dimensional floating point array
  94. ;     containing the interpolated surface, sampled at the grid points.
  95. ;
  96. ; RESTRICTIONS:
  97. ;    The accuracy of this function is limited by the single precision
  98. ;    floating point accuracy of the machine.
  99. ;
  100. ;        SAMPLE EXECUTION TIMES  (measured on a Sun IPX)
  101. ;    # of input points    # of output points    Seconds
  102. ;    16            676            0.19
  103. ;    32            676            0.42
  104. ;    64            676            1.27
  105. ;    128            676            4.84
  106. ;    256            676            24.6
  107. ;    64            256            1.12
  108. ;    64            1024            1.50
  109. ;    64            4096            1.97
  110. ;    64            16384            3.32
  111. ;
  112. ; PROCEDURE:
  113. ;    A minimum curvature spline surface is fitted to the data points
  114. ;    described by X, Y, and Z.  The basis function:
  115. ;        C(x0,x1, y0,y1) = d^2 log(d^k),
  116. ;    where d is the distance between (x0,y0), (x1,y1), is used,
  117. ;    as described by Franke, R., "Smooth interpolation of scattered
  118. ;    data by local thin plate splines," Computers Math With Applic.,
  119. ;    v.8, no. 4, p. 273-281, 1982.  k = 1 for minimum curvature surface,
  120. ;    and 2 for Thin Plate Splines.  For N data points, a system of N+3 
  121. ;    simultaneous equations are solved for the coefficients of the 
  122. ;    surface.  For any interpolation point, the interpolated value
  123. ;    is:
  124. ;      F(x,y) = b(0) + b(1)*x + b(2)*y + Sum(a(i)*C(X(i),x,Y(i),y))
  125. ;
  126. ;    The results obtained the thin plate spline (TPS) or the minimum
  127. ;    curvature surface (MCS) methods are very similar.  The only 
  128. ;    difference is in the basis functions: TPS uses d^2*alog(d^2),
  129. ;    and MCS uses d^2*alog(d), where d is the distance from 
  130. ;    point (x(i),y(i)).
  131. ;
  132. ; EXAMPLES:
  133. ; Example 1: Irregularly gridded cases
  134. ;    Make a random set of points that lie on a gaussian:
  135. ;      n = 15                ;# random points
  136. ;      x = RANDOMU(seed, n)
  137. ;      y = RANDOMU(seed, n)
  138. ;      z = exp(-2 * ((x-.5)^2 + (y-.5)^2))      ;The gaussian
  139. ;
  140. ;     get a 26 by 26 grid over the rectangle bounding x and y:
  141. ;      r = min_curve_surf(z, x, y)    ;Get the surface.
  142. ;
  143. ;     Or: get a surface over the unit square, with spacing of 0.05:
  144. ;      r = min_curve_surf(z, x, y, GS=[0.05, 0.05], BOUNDS=[0,0,1,1])
  145. ;
  146. ;     Or: get a 10 by 10 surface over the rectangle bounding x and y:
  147. ;      r = min_curve_surf(z, x, y, NX=10, NY=10)
  148. ;
  149. ; Example 2: Regularly gridded cases
  150. ;    Make some random data
  151. ;      z = randomu(seed, 5, 6)
  152. ;
  153. ;    interpolate to a 26 x 26 grid:
  154. ;      CONTOUR, min_curve_surf(z, /REGULAR)
  155. ;
  156. ; MODIFICATION HISTORY:
  157. ;    DMS, RSI, March, 1993.  Written.
  158. ;    DMS, RSI, July, 1993.   Added XOUT and YOUT, XPOUT and YPOUT.
  159. ;-
  160.  
  161.  
  162. ON_ERROR, 2
  163. s = size(z)        ;Assume 2D
  164. nx = s[1]
  165. ny = s[2]
  166.  
  167. reg = keyword_set(regular) or (n_params() eq 1)
  168.  
  169. if n_elements(xgrid) eq 2 then begin
  170.     x = findgen(nx) * xgrid[1] + xgrid[0]
  171.     reg = 1
  172. endif else if n_elements(xvalues) gt 0 then begin
  173.     if n_elements(xvalues) ne nx then $
  174.         message,'Xvalues must have '+string(nx)+' elements.'
  175.     x = xvalues
  176.     reg = 1
  177. endif
  178.  
  179. if n_elements(ygrid) eq 2 then begin
  180.     y = findgen(ny) * ygrid[1] + ygrid[0]
  181.     reg = 1
  182. endif else if n_elements(yvalues) gt 0 then begin
  183.     if n_elements(yvalues) ne ny then $
  184.         message,'Yvalues must have '+string(ny)+' elements.'
  185.     y = yvalues
  186.     reg = 1
  187. endif
  188.  
  189. if reg then begin
  190.     if s[0] ne 2 then message,'Z array must be 2D for regular grids'
  191.     if n_elements(x) ne nx then x = findgen(nx)
  192.     if n_elements(y) ne ny then y = findgen(ny)
  193.     x = x # replicate(1., ny)    ;Expand to full arrays.
  194.     y = replicate(1.,nx) # y
  195.     endif
  196.  
  197. n = n_elements(x)
  198. if n ne n_elements(y) or n ne n_elements(z) then $
  199.     message,'x, y, and z must have same number of elements.'
  200.  
  201. m = n + 3            ;# of eqns to solve
  202. a = fltarr(m, m)        ;LHS
  203.  
  204. ; For thin-plate-splines, terms are r^2 log(r^2).  For min curve surf,
  205. ; terms are r^2 log(r).
  206. if keyword_set(tps) then k = 1.0 else k = 0.5
  207.  
  208. for i=0, n-2 do for j=i+1,n-1 do begin
  209.     d = (x[i]-x[j])^2 + (y[i]-y[j])^2 > 1.0e-8  ;Distance squared...
  210.     d = d * alog(d)* k    ;d^2 * alog(d^(1./k))
  211.     a[i+3,j] = d
  212.     a[j+3,i] = d
  213.     endfor
  214.  
  215. a[0,0:n-1] = 1        ; fill rest of array
  216. a[1,0:n-1] = reform(x,1,n)
  217. a[2,0:n-1] = reform(y,1,n)
  218.  
  219. a[3:m-1,n] = 1.
  220. a[3,n+1] = reform(x, n, 1)
  221. a[3,n+2] = reform(y, n, 1)
  222.  
  223. b = fltarr(m)        ;Right hand side
  224. b[0] = reform(z,n,1)
  225.  
  226. c = b # invert(a)    ;Solution using inverse
  227. ;  c = svd_solve(transpose(a),b)    ;solution using svd
  228.  
  229. if n_elements(XPOUT) gt 0 then begin    ;Explicit output locations?
  230.   if n_elements(YPOUT) ne n_elements(XPOUT) then $
  231.     message, 'XPOUT and YPOUT must have same number of points'
  232. endif else begin            ;Regular grid
  233.   if n_elements(bounds) lt 4 then begin    ;Bounds specified?
  234.     xmin = min(x, max = xmax)
  235.     ymin = min(y, max = ymax)
  236.     bounds = [xmin, ymin, xmax, ymax]
  237.     endif
  238.  
  239.   if n_elements(gs) lt 2 then begin    ;GS specified?  No.
  240.     if n_elements(nx0) le 0 then nx = 26 else nx = nx0 ;Defaults for nx and ny
  241.     if n_elements(ny0) le 0 then ny = 26 else ny = ny0
  242.     gs = [(bounds[2]-bounds[0])/(nx-1.), $
  243.        (bounds[3]-bounds[1])/(ny-1.)]
  244.   endif else begin            ;GS is specified?
  245.     if n_elements(nx0) le 0 then $
  246.     nx = ceil((bounds[2]-bounds[0])/gs[0]) + 1 $
  247.     else nx = nx0
  248.     if n_elements(ny0) le 0 then $
  249.     ny = ceil((bounds[3]-bounds[1])/gs[1]) + 1 $
  250.     else ny = ny0
  251.   endelse
  252.  
  253.  
  254.   if n_elements(xout) gt 0 then begin        ;Output grid specified?
  255.     nx = n_elements(xout)
  256.     xpout = xout
  257.   endif else xpout = gs[0] * findgen(nx) + bounds[0]
  258.  
  259.   if n_elements(yout) gt 0 then begin
  260.     ny = n_elements(yout)
  261.     ypout = yout
  262.   endif else ypout = gs[1] * findgen(ny) + bounds[1]
  263.   xpout = xpout # replicate(1.,ny)
  264.   ypout = replicate(1., nx) # ypout
  265. endelse                    ;Regular grid
  266.  
  267. ; For min_curve_surf, divide c[3:*] by 2 cause we use d^2 rather than d.
  268. if keyword_set(tps) eq 0 then c[3] = c[3:*]/2.0
  269.             ;common scale factor
  270. s = c[0] + c[1] * xpout + c[2] * ypout        ;First terms
  271. for i=0, n-1 do begin
  272.     d = (xpout-x[i])^2 + (ypout-y[i])^2 > 1.0e-8  ;Distance squared
  273.     s = s + d * alog(d)* c[i+3]
  274.     endfor
  275. return, s
  276. end
  277.